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Abstract 

The one-dimensional forest-fire model including lightnings is studied numerically and ana- 
lytically. For the tree correlation function, a new correlation length with critical exponent 
z/ ^ 5/6 is found by simulations. A Hamiltonian formulation is introduced which enables 
one to study the stationary state close to the critical point using quantum-mechanical per- 
turbation theory. With this formulation also the structure of the low-lying relaxation spec- 
trum and the critical behaviour of the smallest complex gap are investigated numerically. 
Finally, it is shown that critical correlation functions can be obtained from a simplified 
model involving only the total number of trees although such simplified models are unable 
to reproduce the correct off-critical behaviour. 
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1. Introduction 



Power laws and scaling behaviour, familiar from equilibrium critical phenomena, can also 
be found in non-equilibrium systems. If in such a case no fine-tuning of a parameter is 
necessary, the situation has been described as self-organized criticality [1,2]. To illustrate 
the phenomenon, simple models for sandpiles [1,2], forest fires [3] or evolution [4,5] have 
been proposed. Their common feature are avalanche-type processes in the dynamics. In 
the case of the forest- fire models (fi'm) , the originally proposed version did not show proper 
scaling behaviour [6, 7], and it was necessary to introduce lightning strokes in order to find 
it [8]. Even then, the spatial power laws are valid only up to a typical length depending 
on the lightning rate. Therefore the ffm is in general not at a critical point, but only close 
to it. This is also the picture obtained from a renormalization treatment [9]. Nevertheless, 
the ffms are interesting systems (for a brief review see [10] and for a recent account [11]), 
and a detailed understanding of their properties, in particular in comparison with usual 
critical systems, is desirable. This should be simplest in the one-dimensional case. 

In one dimension, a few exact results have been obtained [12,13]. In partiuclar, the 
stationary distribution of tree clusters with size s is given by [13] 

n(s) = , , (1.1) 

^ ' (s + l)(s + 2) ^ ' 

where p is the density. This expression is valid for s <^ and gives n[s) ~ for large 
s. The length depends on the rates p (for tree growth) and / (for lightning strokes) 
as ~ vl f 1 up to logarithmic corrections. Thus, the corresponding exponent \s> v = \. 
The time correlation function of burning trees could also be obtained. These results were 
checked by numerical calculations and a certain overall picture has emerged. However, the 
precise nature of the stationary state or the form of the relaxation spectrum are still to be 
determined. 

To make some progress in this direction, we first studied the stationary correlations in 
more detail. In particular, we looked at the simple tree correlation function and found an 
unexpected result, namely a different correlation length which diverges with an exponent 
vt ^ 5/6. This is in contrast to results in two dimensions [14] and indicates a rather 
complicated structure of the stationary state. To obtain this state explicitly, one has to 
treat the master equation. For this we used a quantum-mechanical formulation which 
has proved to be quite useful in other problems (see for example [15]). If one assumes 
instant burning of tree clusters, one then obtains a spin one-half quantum chain with both 
single and cluster-fiip processes. A number of eigenstates can be found exactly, but in 
general one has to resort to numerics. We calculated the low- lying spectrum for systems 
up to L = 20 sites, obtaining various branches in the complex plane and a gap which closes 
algebraically as the rate / goes to zero. We also considered the limit / ^ analytically and 
found that to order 0(/) the ground state lies in a subspace of L -|- 1 totally symmetric 
configurations. Working only in this subspace, one obtains a simple model where the 
system evolves in cycles [16]. Then the stationary state can be found explicitly and gives 
the cluster distribution (1.1) but it contains no correlations. Also the low-lying spectrum 
differs from the exact one. A similar observation for models of evolution has already been 
made in [17] where it was concluded that the behaviour of only one quantity is not a 
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sufficient criterion for criticality. To pursue this aspect further, we also investigated to 
what extent the cluster-size distribution of the full model for / > can be recovered in 
the symmetric subspace. 

2. Simulations of correlation functions 

This section presents results of simulations of two correlation functions. First, the two- 
point function of trees is studied and it is shown that the exponent of the corresponding 
correlation length is not an integer. By contrast, the correlation function inside a single 
cluster shows a critical behaviour which is consistent with the critical exponent of the 
cluster-size distribution n{s) eq. (1.1). 

The ffm introduced in [3] is defined on a cubic lattice in d dimensions. Any site can have 
three states: It can be empty, or it can either be occupied by a green tree or a by a burning 
tree. The dynamics of the model was specified in [3] by the following parallel update rules 
from one timestep to the next: a) A burning tree becomes an empty site, b) A green 
tree becomes a burninng tree if at least one of its nearest neighbours is burning, c) At 
an empty site a green tree grows with probability p. In order to obtain proper critical 
behaviour, the following rule was added in [8]: d) Every green tree can spontaneously be 
struck by lightning with probability / and thus become a burning tree even if no neighbour 
is burning. 

Quantities of interest in this model are in particular the forest clusters, i.e. maximal con- 
nected sets of green trees. In one dimension a forest cluster is a string of green trees 
bounded by empty sites and its size s is just the number of trees. 

The critical point of the above ffm, which has been studied as an example of self-organized 
criticality, arises at p — > and //p — > 0. Keeping f/p fixed and taking p — > has the effect 
that once a cluster is struck by lightning it burns down before anything else happens. After 
redefining the time scale, the problem reduces to a two-state model where only empty sites 
(E) and trees (T) occur. The dynamics of this effective model is specified by the following 
parallel update rules: 

1) At an empty site a tree grows with probability p. 

2) A forest cluster of size s is struck by lightning with rate fs and all trees inside it 
become empty sites during one timestep. 

Simulations have been carried out almost always with this reduced model (see e.g. [18, 14]). 
Below we will always think of this reduced version if we refer to 'the forest-fire model'. 
We will restrict our attention to the one-dimensional version with periodic boundary con- 
ditions. The size of the lattice will be denoted by L. 

First we have investigated two-point functions of trees (T^Tx+y) in one dimension. To this 
end we have performed simulations for 0.0025 < p < 0.01 in the range 2 10~^ f /p 
2 10~^. The sizes of lattices ranged between L = 4000 for large / and L = 500000 for small 
/. One random initial condition with density p = | was chosen. Then to iterations were 
performed in order to equilibrate the system (usually to = 3000, only for very small / one 
needs a few more iterations). After that {TxT^+y) was determined every 100 timesteps 
from to to to + 1900 by averaging over all x. This determination of (T^Tx+y) is much more 
time consuming than the simulation itself. The procedure was repeated 10 times with the 
same initial conditions, amounting to a total of 200L measurements for {TxT^+y). 
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The resulting data can nicely be fitted by 



C{y) 



(TxTx+y) 



-{Tx 



ae 



-y/i 



(2.1) 



for 1 < y. (Note that (T^) ~ p and = T^.) The parameters a and ^ were estimated 
by taking the logarithm of the r.h.s. of (2.1) and then performing a linear regression for 
1 < y < I/max- ymax was chosen such that statistical errors can be neglected for y < j/max- 
With the statistics of our simulations j/max was always a few correlation lengths. 
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Fig. 1: The correlation function C{y) for L = 100000, p = 1/200, 
f /p = 1/100 (full line). The dashed line is the exponential form (2.1) 
with a = 0.0707, C = 31.81. 



Fig. 1 shows the result of a detailed simulation where the average was taken every 100 
timesteps over the larger time interval from to = 3000 to to + 39900. The agreement of 
the simulation with the exponential form is very good for a few correlation lengths. For 
y ^ 200 statistical errors start to become so large that the difference on the l.h.s. of (2.1) 
cannot be determined sufficiently accurately any more. For smaller values of y there are 
small wiggles over longer ranges of the distance y which are typical for the simulations and 
seem to be caused by statistical fluctuations. Apart from that, no substantial corrections 
to the form (2.1) are visible. 

Fig. 2 shows the results of the correlation length ^ obtained from simulations as a function 
of f /p. There are some residual statistical as well as systematic errors due to a small p- 
dependence which are, however, not larger than the symbols. One can see that the values 
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for ^ (and similarly for a) are in good agreement with the form 



a ~ 



(2.2) 



Performing linear regression fits on a doubly logarithmic scale one finds 



i^T = 0.8336 ± 0.0036 , /i = 0.1031 ± 0.0022 . 



(2.3) 



Surprisingly ut agrees very well with the gap index z/ = 5/6 for the three-state Potts model 
[19]. The most probable rational value for the other exponent is n — 1/10. 
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Fig. 2: The correlation length ^ as a function of f/p for various values 
of /, p and L. 



It was argued in [14] for two dimensions that the critical exponents are independent of 
the particular correlation function considered. In order to investigate this question in one 
dimension we also looked at the correlation function {TxTx^y)c describing the probability 
to find two trees at positions x and x + y inside the same cluster. In d = 1 one has 

) =■■ K{y) . (2.4) 

This quantity can easily be simulated as before. Its determination is actually much less 
time consuming. 
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The correlation function K{y) can be fitted with a single exponential function for y suffi- 
ciently large, but for smaller y an exponential function is not a good approximation. Thus, 
we fit K{y) by 

oo 

K{y) = J2 a^^^e-^'/^i'^^ (2.5) 

where the ^jf^ decrease with r. The first few terms in the sum (2.5) already give excellent 
approximations. In an interval T/min < ?/ < ?/max all terms with r > 1 can be neglected and 
this can be used to estimate the largest decay length = Cc^^ with a linear regression for 
the logarithm of the correlation function K{y). Typically one finds ymin ~ Cc, J/max > 5^c- 

Note that K{y) is related to the cluster-size distribution n{s) via 

s>y 

This can be inverted to give 

n{s) = K(s-1) -2K{s) + K{s + l). (2.7) 

This implies that if K{s) decays exponentially in some region, also n{s) must decay expo- 
nentially in roughly the same region. However, the lattice Laplacian in (2.7) increases the 
weight of the next to leading terms in (2.5). Therefore, the purely exponential decay is ex- 
pected to set in at larger s for n{s) than for K{s). The value sq where a single exponential 
function starts to describe n{s) can be estimated as follows. One requires that the form 
n{s) — (1 — p)/s^ can be matched up to the linear term at sq with a single exponential func- 
tion for K{s), i.e. K{s) — ae~^l^'' . According to (2.7) one must have n(so) = ^-ftr(s)|s=s„ 

and ^?^(s)|s=so = ■^-^('5)|s=S()- This leads to the relation a ^ {\ — p) (which is indeed 
well verified in simulations) for the normalization constants and sq ~ the desired 

crossover point. 

Fig. 3 shows the correlation function i^(s) and the normalized cluster-size distribution 
n{s). The simulation was performed in the same manner as in Fig. 1. n{s) was determined 
by counting all clusters of size s in the system during each of the 40000 timesteps after 
to = 3000. In this time-interval, estimates for K{\j) were obtained every hundreth timestep 
as a spatial average over the entire system. The huge amount of data in particular for n(s) is 
needed in order to obtain small statistical fiuctuations without too much extra smoothing. 
One observes that both expectation values decay exponentially over a large range of s. 
The correlation lengths obtained by fits in this region agree very well with each other as 
it should be according to (2.7). The actual value is ~ 28. Furthermore, the crossover of 
the cluster-size distribution from the form (1.1) to the exponential decay does indeed occur 
in the vicinity of one correlation length, the exponential decay being well visible beyond 
s = lO^c- It seems that K{y) starts to decay faster than exponentially around s = 400 in 
Fig. 3. However, comparison with simulations with a smaller amount of samples indicates 
that this effect is a residual statistical error. Thus, the exponential decay of n{s) and K{s) 
could well extend indefinitely, the probability of clusters larger than 450 trees simply being 
so small that they did not appear in the present simulation. 



5 



Although the simulation in Fig. 3 was performed for a fairly large value of //p, it does 
reflect the typical situation also closer to the critical point f /p = 0. The form (1.1) for n{s) 
is valid approximately up to one correlation length while the exponential decay can be 
verified for at least a few ^c- For smaller f /p than in Fig. 3 there is a small difference: There 
is a region where the form (1.1) for n{s) lies below the actual exponential decrease. This 
corresponds to the 'bump' observed in a logarithmic plot of the cluster size distribution 
[13]. 
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Fig. 3: The correlation function K{s) and the cluster-size distribution 
n{s) for L = 100000, p = 1/200, f/p = 1/100. 

There is a striking similarity of n{s) in Fig. 3 with the avalanche distribution shown in 
Fig. 6b of [20] for a different two-dimensional forest-fire model. This similarity might 
indicate that one could expect similar results for the two-dimensional version of the forest 
fire model considered in this paper, but it may also just illustrate the fact that cluster size 
distributions of very different models can qualitatively resemble each other. 

We have also performed simulations to study the //p-dependence of K{y). They were 
done for 4 10"^ < f/p <2 10"^ and 15000 <L< 500000. One finds that the resulting 
data is consistent with 

In(ec) = (0.95 ± 0.06) ^ . (2.8) 

The functional form of the relation (2.8) was already given in Ref. [13] although the 
derivation has to be considered with some care (see section 5 below). 

This result is to be contrasted with the one for the two-point correlation function C{y). 
The two-point function is significantly more accurately descibed by a single exponential 
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function (compare Figs. 1 and 3). Even more, the critical exponents are clearly different 
- see eq. (2.3). Thus, in one dimension there seems to be no direct relation between the 
two different correlation functions. Note that the relative errors in Ref. [f 4] for the critical 
exponents in two dimensions are much larger than the ones of our simulations. Therefore it 
is conceivable that a more precise study could reveal different exponents in two dimensions 
as well. 

Finally, we would like to mention than an exponential decay of n(s) as found above, 
can already be obtained from a mean-field approximation. Complete factorization of the 
correlator defining n{s) yields n{s) = (1 — p^p^, i.e. n{s) varies exponentially. In the 
mean-field approximation for the original version of the ffm with burning trees one finds 
the critical behaviour 1 — p ~ f /p. This implies = ~l/lii(p) ^ p/ fi i-e. the critical 
exponent is = 1 in mean-field approximation. This agrees with (2.8) up to logarithmic 
corrections but one should note that the critical behaviour of p is not right. The correct 
result is p/ (f — p) ~ ln(p//) [13] and the corresponding exponent differs by one from the 
mean- field result. 

3. Hamiltonian formulation of the model 

In this section we present a quantum-mechanical formulation of the ffm and use it to 
discuss correlations at the critical point. 

In order to set up notations, let us briefly recall the Hamiltonian formulation of stochastic 
lattice models (see e.g. [15]). Denote the probability to find the system in a configuration 
by P({/3}). Then the time evolution of the system is given by the master equation 

dtPm) = - E W{p}^{p,^pm) + E W{p'}^{p^p{m) (3.1) 

where the two terms on the r.h.s. represent the loss and gain processes respectively. The 
symbol dt can equally well denote evolution in discrete or continuous time. To formulate 
(3.1) in quantum-mechanical terms one introduces a basis of states | {/?}) corresponding 
to the configurations {/?} and specifies the state of the system by a vector 

m 

Then the master equation (3.1) takes the form of a Schrodinger equation in imaginary time 

dt\P)=-H\P) (3.3) 

where the Hamiltonian H for a stochastic system is in general non-hermitean. Due to 
probability conservation one always has an eigenstate of H with eigenvalue zero. 

In order to apply this general formalism to the ffm we choose the following basis at each 
site: 

^ ^ ^ = empty place = , ) ~ ^^^^ ~ ^ ' i'^-^) 
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A general state | {/?}) corrsponding to a configuration of the complete system is given by 

\m) = |^l...^L) , (3.5) 

where = 0,1 encodes the state of the site x. The vectors (3.5) form a basis of a 2^ 
dimensional complex vector space. The notation and Ej; will be used for the operators 
that measure if site x is occupied by a tree or empty. 

We will consider the case of small p and /. Then in each time step changes occur essentially 
only at one site and one can replace parallel dynamics by sequential dynamics. Note that 
for the original three-state model parallel dynamics is essential because e.g. the probability 
of fire spreading simultaneously at different places is large. With sequential dynamics the 
Hamiltonian takes the form 

H = pHo + fV (3.6a) 

and can be written down explicitly using Pauli matrices. In terms of Pauli matrices the 
operators and are given by T-j; = ^ (1 + a^)^ and E^ — ^ {1 — cr^)^. The state of a 
given site is changed by the operators c = | (c^ ± ia^). 

Hq is a sum of single-site terms (i^o)x describing growth of a tree at site x. In spin 
language, the gain term then involves single spin flips from to 1: 



^0 = ^(^1-^^).-^.+) • (3-6b) 

x=l 

The term V describes the burning down of tree clusters and is given by 



L L 



^ = E E M 1(1 - \ (1 + •••1(1+ cr^Ui 1(1- cr^Ui^, 

-1(1 - a^)^ a-+, . ■ ■ a-^, \ (1 - a^),+,+ J 



(3.6c) 



with suitable conventions such that the projectors at the boundary drop out for a cluster 
that occupies the entire system. This term contains flips of entire clusters from state 1 to 
the empty state 0. The weight of such a cluster-flip process is proportional to the size I of 
the cluster. 

We would like to note that our Hamiltonian (3.6) is different from the one presented in 
Anhang B of [21]. Also the Fock space formulations of the master equation in [22,23] are 
different since they treat situations with more than two states per site. 

Clearly, it is possible to choose p = 1 in (3.6a) by suitably rescaling H. Below we will 
assume p = 1 and we will be interested in small /. Then one can apply standard quantum- 
mechanical perturbation theory to (3.6). 

Before performing any computation one can already see that if one wants to take the limit 
L — > oo of a perturbation expansion, one is forced to approach the critical point / = 0. 
This is due to the fact that a perturbation series is expected to converge as long as / times 
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the largest matrix element of V (which has value L) is small compared to the distance 1 
between the eigenvalues of Hq. This amounts to imposing 



/ L « 1 (3.7) 

in order to be on safe grounds. Thus, one is forced to let / — > as L — > oo. In other words, 
the radius of convergence of the perturbation series is expected to be zero for L = cx). Still, 
perturbation theory will tell us something about the properties of the critical point itself. 

First, we use the Hamiltonian formulation of the ffm in order to determine the stationary 
state perturbatively. In order to simplify the presentation we introduce the notation 

\N):^Y.\ ^ ) (3.8) 

TV empty places 

for a sum over all configurations with a total of empty places. For the stationary state 
I G) perturbation theory is particularly simple. One makes the power-series ansatz | G) = 
J2T=o \'^'^) ■ Obviously |Go) = 1 1 • • • 1) = |0). For the first-order correction jCi), 
the standard quantum- mechanical formulas lead to the condition Hq \Gi) +V\Go) — 
where the projection onto | Go) has to be subtracted. Since Hq as well as V \ Gq) are 
invariant under permutations of the spins, \Gi) must also have this property and thus 
can be written as a superposition of the states | N) 

L 

\Gi) = J2^n\N) . (3.9) 
Ar=i 

Let us note that, apart form the growth processes, only one lightning process enters in 
this first-order calculation, namely the one leading from the completely full system to a 
completely empty system. 

Using now that the coefficient of a state | N) with 1<A'"<L — lin the expression for 
Ho\Gi) must be zero, one finds the recurrence relation qjtv+i = N/ (L — N) aN ior 1 < N. 

The only further condition is ckl = 1. This gives ckjv = (^Zi) • In summary, we have 
found that 

\G) = |1. ..!)+/ ETZ^|iV)+W) (3.10) 
statisfies I G) =0(f). 

Eq. (3.10) is in the form that arises naturally from perturbation theory, but is is not yet 
properly normalized as a probablility. Thus, we multiply the state | G) with A/l and 
impose the normalization condition that the sum of all coefficients in the state A/l | G) is 
one. This leads to 

^=1+/ i:Ms+o{f)=i+fLf2j^+o{f). (3.11) 

N=l \N-i) 7V=1 
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From (3.10) one can compute every correlation function for / — > 0. Let 

riris) = {E,^ . . • • (3.12) 

be any correlation function with s trees and r + 1 empty places at certain fixed sites. Since 
(3.10) is invariant under a permutation of the L sites of the lattice, this correlation function 
will only depend on the total number of trees and empty places, respectively. Prom (3.10) 
one finds that 

for any finite system, r > and 0<s<L — r — 1. The second identity is a purely 
combinatorial one and proven in appendix A for r — 0. For r > the result follows most 
easily from the equation of motion for nr{s) in the stationary state (compare section 5). 
Note that any correlation function involving at least one empty place decays algebraically 
in this approximation, the chain length L acting as cutoff. Furthermore, tree-tree or 
cluster-cluster correlation functions are trivial in the sense that they do not depend on the 
distance between the trees, respective clusters. 

Let us briefly comment on a few special cases of (3.13): 

1) For the cluster-size distribution n{s) — ni{s) one recovers the result eq. (1.1) with 
1 — p = MlJL. a derivation of this result for r = 1 along similar lines as above and in 
section 5 below can also be found in [13, 16]. In fact, the discussion of the stationary 
state at small / in [13] inspired us to apply quantum-mechanical perturbation theory 
to the ffm. 

2) The probability to find a cluster of s empty places with two adjacent trees is given by 
ns_i(2) = iMhf^l {s{s -I- l)(s -I- 2)). The functional form of this correlation function 
was already derived in [12]. Note that in eq. (5) of [12] which was used to obtain this 
result, the loss terms by a lightning stroke at the boundaries are missing. This is, 
however, irrelevant for the derivation of the form of ?is_i(2). 

3) For the probability to find a string of s adjacent empty sites one finds that ns-i(O) = 
MhfL/s for 1 < s < L — 1. For this special case it is also comparably easy to directly 
prove the combinatorial identity (3.13) by induction on L. 

4) The probability K{s) to find s trees at certain places and no empty places can be 
obtained from (2.7) using (3.13) for n{s). One finds in particular that K{s) ~ ln(s) 
for large s. Note that the smaller the exponent of a correlation function the sooner 
will it be dominated by an exponential decay for f ^ 0. In particular, the asymptotic 
logarithmic behaviour of K{s) is therefore not observable in simulations (compare Fig. 
3 of section 2). 

It is also straightforward to compute the mean density of trees p. One finds that p/ (1 — p) ~ 
ln(L) for L large. This also shows that in this first-order approximation only the finite size 
of the system prevents it from being critical. 

We now proceed with a discussion of the relaxation spectrum. For / = the operator 
H — Hq has real eigenvalues which are simple integers A = 0, 1, . . . , L. In general, one 
can use the fact that for periodic boundary conditions the lattice translation operator 
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commutes with the Hamiltonian and the spectrum sphts into L sectors with different 
momenta P. We introduce the notation 




for states of fixed momentum P = 27rn/L and specify the first Brillouin zone by < n < L, 
n e Z. 

Some of the excited states can easily be written down. Using that a sum of all Lth roots 
of unity vanishes, one can see that 



for P 7^ 0. Eq. (3.15) implies that there are exact excited states with eigenvalue 1 + /(L — 1) 
and 2 + /(L — 2) respectively. However, for L large they will not belong to the low-lying 
part of the excitation spectrum. 

Other excited states can be treated with standard quantum mechanical perturbation the- 
ory, bearing in mind that the left and right eigenvectors have to be treated separately. 
Some computations are presented in appendix B. The final result for the states with eigen- 
value A = 1 of i^o is that H has L — 1 times the eigenvalue 1 -|- /(L — 1) + 0{p) and 
one times the eigenvalue 1 + /(2L — 1) + 0{f'^). Similarly, the states with eigenvalue 
A = 2 of Hq give rise to an eigenvalue 2 + /(L — 2) + 0{f'^) of H for each momentum 
(including the translationally invariant sector). In the translationally invariant sector this 
eigenvalue 2 + f{L-2)+0{f'^) can cross with the eigenvalue 1 + f{2L - 1) + 0{p) around 
/ = 1/(L-|-1). The physically interesting region is located beyond this crossing where 
the lowest eigenvalues form a complex-conjugate pair. However, there are fundamental 
obstacles to reach this region beyond the crossing with perturbation expansions around 



4. Numerical computation of the relaxation spectrum 

In order to obtain results for the thermodynamic limit of the relaxation spectrum, one has 
to resort to a numerical diagonalization the Hamiltonian at finite L. In this section we 
argue on the basis of such computations that the ffm has a non-vanishing complex gap in 
the thermodynamic limit for / > 0. The exponents describing the critical behaviour for 
/ — s> of the real and imaginary parts of this complex gap are both non-integral. 

After projecting onto the eigenspaces (3.14) with fixed momentum P, the Hamiltonian 
(3.6) becomes a square matrix of approximate size 2^/L. Up to L = 20, where the size of 
the matrix is slightly larger than 50000, a few extremal eigenvalues can be obtained using 
standard methods. 

In addition to translational invariance, the model is invariant under parity x i— > —x. This 
can be exploited to reduce the dimensionality of the problem a little further for P = and 
P = TT. More important are the following consequences of the invariance with respect to 
parity: 

1) The spectra in the sectors with momenta P and 2tt — P are equal. 



i?||01...1))p=(l + /(L 
iy||001...1)) p =(2 + /(L 



1) )||01...1))p 

2) )||001...1))p + (l + e^O||01...1))p 



(3.15) 



/ = o. 
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2) In each sector with fixed momentum, the eigenvalues are either real or come in complex 
conjugate pairs. This holds because the characteristic polynomial of the complete 
Hamiltonian (3.6) is real, and because of property 1). 

Fig. 4 shows the low-lying excitations A for L = 20 and / = 0.1 in the complex plane. 
Because of the symmetries mentioned above we have restricted ourselves to < P < tt 
and to Im(A) > 0. Note that the figure does not include the ground state A = and that 
there are many more eigenvalues near its right end and beyond. 
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Fig. 4: Eigenvalues of the forest-fire Hamiltonian with L = 20 sites, 
p = 1, / = 0.1. For the symbols compare the text. 



The comparably large value of / = 0.1 has been chosen in order to have the dispersionless 
modes (3.15) separated from the low-lying part of the spectrum. The first one of these 
modes is located at A = 2.9 in Fig. 4. One observes a group of levels on the real axis 
indicated by the symbols '-|-'. The level with the smallest gap has P = Att/L. Then the 
gap increases with increasing P and reaches its maximum at P = tt. One can also clearly 
see another sequence of complex- conjugate levels (which are drawn with the symbol 'O') 
starting around A = 1.4 ± 2i for P = and tending to A = 2.6 ± 0.56i for P ^ tt. 
There are also further levels in this region. The ones we have computed are marked by 
the symbol '□' in Fig. 4. They do not form groups as clearly as the levels described 
before. It is also possible that some of them will still move considerably with increasing 
L due to their interaction with the mode at A = 1 4- f{L — 1). The figure demonstrates 
that the relaxational spectrum of the Hamiltonian (3.6) does not have a simple particle 
interpretation. The levels drawn with the symbols and 'O' are isolated levels for P fixed 
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and therefore clearly not multi-particle states which would form continua, neither are their 
features typical for single-particle states. In particular, single-particle states usually have 
exponentially small corrections in the chain length, but the levels under consideration have 
larger corrections (compare also below). Even if the levels corresponding to the symbols 
'-h' and '<>' were considered as fundamental single-particle states, some of the remaining 
levels which we denoted by '□' would remain unexplained. 

In the remainder of this section we concentrate on the behaviour of the lowest gap Agx in 
the translationally invariant sector. In Fig. 4 one has Aex ~ 1.3579 ± 1.9586z. Whether 
the level Aex actually has the smallest real part, or if other levels exhibit different critical 
behaviour would require more detailed investigations which are beyond the scope of the 
present paper. The restriction to P = is convenient because only in the translationally 
invariant sector the finite-size effects can be analyzed easily. In addition, measurements 
involving spatial averages project onto this sector, so that it plays a particular role. 



Table 1 contains asymptotic values obtained numerically for the gap Agx- They were 
obtained by computing numerically the lowest complex conjugate levels ^) for 5 < L < 20 
and extrapolating to L = oo using the van den Broeck-Schwartz algorithm (see e.g. [24]). 





diagonalization {j> — 1) 


simulation with 


p = 1/200, L = 10000 


f/p 


lim Aex 


error of lim Aex 


relaxation time 


oscillation period 




L— »oo 




0.1 


1.2947 ±2.0030z 


0.0003 


140 


560 


0.09 


1.2344 ± 1.9625z 


0.0003 


150 


575 


0.08 


1.1726 ±1.9170z 


0.0002 


160 


615 


0.07 


1.1069 ±1.8664z 


0.0002 


160 


695 


0.06 


1.0360 ± 1.8093z 


0.0004 


190 


705 


0.05 


0.9600 ±1.7461z 


0.0004 


205 


715 


0.04 


0.8767 ±1.6721z 


0.0003 


255 


740 


0.03 


0.7820 ± 1.5820i 


0.0010 


270 


775 


0.02 


0.6780 ±1.4680i 


0.0040 


300 


875 



Table 1: Estimates for the first translationally invariant gap of the forest-fire Hamiltonian. 



The values in table 1 have to be considered with some care because the maximal system 
size L = 20 used for the extrapolation to L = oo is still quite small. However, all sequences 
used for the extrapolation are monotonuos in L. Furthermore, the deviation at finite L 
from the values given above is typically of the order L~^. Under these conditions one 
can usually obtain reliable estimates for the thermodynamic limit L — oo although the 
estimates in table 1 for the error of Aex may be somewhat too optimistic. 

In order to perform a further check of the extrapolations we have simulated the temporal 
behaviour of the average density of trees p{t). This quantity approaches its limiting value 
p(oo) in an oscillatory way. The simulations were performed for p = 1/200 on a lattice of 



^) Typically these levels fall below others (with respect to the real part), which are 
located on the real axis, for L ~ 10 in the range of / considered here. 
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size L = 10000 - a size that is much closer to the thermodynamic hmit. For one random 
initial condition with p(0) = | the time evolution of p{t) was averaged over 100 realizations. 
Then the oscillation period was determined by looking for crossings of p{t) at early times 
t with the stationary value p{oo) (which can be estimated from p(t) for large t). The 
relaxation time can be obtained by looking at the values of local extrema of p{t) at early 
times. If only one complex pair contributes, the relaxation time equals l/{p Re(Aex)) while 
the oscillation period is given by '2n/{p Im(Aex)), where Aex is the normalized relaxational 
mode at p = 1. One obtains good agreement with the extrapolations in Table 1 keeping 
the crude method in mind which can at most be expected to be accurate to 10%. 

Because the extrapolations can be performed much more accurately we used them in order 
to determine the behaviour of Agx as a function of /. The data in Table 1 can nicely be 
fitted by 

lim Aex = 3.258 f-^^^ ± 3.127z /"'^^^ . (4.1) 

I/— >oo 

This is consistent with a standard critical point at / = - like the simulations of the 
correlation functions in the previous section. We would like to point out that this critical 
behaviour arises because the limits / ^ and L ^ oo do not commute, Taking f —> 
first, one has Aex = 1 for all L. 

The critical exponents in (4.1) may still have systematic errors due to the fact that only 
systems of maximal size L = 20 were used and that because of this restriction in the 
system size we could not get closer to the critical point than / = 0.02. Thus, the critical 
exponents for the oscillation frequency Im(Aex) and inverse relaxation time Re(Aex) could 
be 1/5 and 2/5, respectively. 

5. Simplified models in the symmetric space 

In this section we consider several simplified models that deal only with the total number 
of trees respectively empty places. The main motivation is that the critical correlation 
functions (3.13) follow most easily from such a simplified model. Because this simplified 
model exhibits neither the critical behaviour (2.8) nor the gap (4.1) we investigate whether 
it is possible to generalize the model in the symmetric space such that it describes also 
some aspects off the critical point correctly. It turns out that the natural generalizations do 
not have these properties. This shows that most of the critical exponents of the ffm seem 
to be due to actual correlations while the power laws of the critical correlation functions 
(3.13) are a global effect. 

The first simplified model corresponds to the first-order expansion (3.10). It consists in 
modifying the dynamics in such a way that this first-order expansion of the ground state 
becomes exact. Since only the burning of the completely full system contributes in first 
order perturbation theory, this is achieved by dropping all lightning processes except for 
the completely full system (see also [16]). This modification leaves the space of states (3.8) 
invariant, so one can restrict oneself to this subspace. We call it symmetric since the states 
are invariant under permutations of the lattice sites. Consequently, there is no real spatial 
structure for this simplified model any more. 

It is useful to replace the basis (3.8) by the following differently normalized one: 

1^):=^- (5.1) 
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The normalization factor in (5.1) is precisely the number of configurations contributing to 
I A^) . In this basis the simplified Hamiltonian is given by the following (L + 1) x (L + 1) 
matrix: 

f fL -1 ••• \ 







1 







(5.2) 



-L 
\-fL ••• Lj 

The normalization of the states in (5.1) is chosen such that the Hamiltonian (5.2) is a 
manifestly stochastic one, i.e. the columns sum to zero. It describes a single particle 
hopping on a ring with sites 0, . . ., L with the following rates: 

(5.3) 



L- 1 



By construction, the state (3.10) is the exact ground state of the simplified Hamiltonian 
(5.2). In the basis (5.1) it reads 



\G) = |0)+/L Eiv'^^- 



(5.4) 



iV=l 



One can use this simplified model to obtain the result (3.13). The equation of motion for 
the quantity nr{s) (see (3.12)) is: 



dfUris) — — (r + l)nr{s) + snr+i{s — 1) 



(5.5) 



for < s < L, r > 0. The first term describes the loss term by tree growth at any of the 
r + 1 empty places whereas the second term is the gain term by a tree having grown at any 
of the s places that are occupied by trees in the final configuration. Note that lightnings 
affect only configurations with s = or s = L trees and are therefore not relevant for the 
present consideration. In the stationary state, eq. (5.5) leads to the recurrence relation 
nr+i{s — 1) = (r + l)nr{s)/s. Together with the combinatorial proof of (3.13) for r = 
presented in appendix A this proves (3.13) for general r. 

Now we proceed with a discussion of the excitation spectrum of Hg. The characteristic 
polynomial of yields the equation 



{fL-X)ll{r-X)-fLllr = 



(5.6) 



r=l 



for the eigenvalues A. A discussion of this characteristic polynomial leads to the following 
results on the relaxational spectrum 



Im(Afc) = 



2nk 



keZ, 



Re(Ai) 



31n(L)^ 



(5.7) 
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for L large. A numerical computation of the lowest eigenvalues of up to L « 10^ yields 

good agreement with the result (5.7) for Im(Afc). The precise behaviour of the real part is 
more difficult to verify, but from the numerical computations one sees that the real parts 
of at least the lowest eigenvalues tend much faster to zero than their imaginary parts. 

The result for Im(A), as well as the eigenfunctions can also be obtained from a continuum 
approximation. Details of the computation are presented in appendix C. The continuum 
model shows that the eigenfunctions are essentially plane waves in the variable ln(A'") / ln(L) 
multiplied by 1/N. In this variable, the kih eigenfunction has k periods. Therefore, the 
value l/Im(AA;) can be interpreted as the average time that a particle needs to hop around 
the /cth part (on this logarithmic scale) of the chain. 

In summary, the gap of the simplified Hamiltonian (5.2) vanishes as l/ln(L) as L — cxo 
irrespective of the value of / ^). This is different from the behaviour of the full ffm as 
observed in the previous section. 

The simplified model (5.2) discussed so far is not able to describe the behaviour of the full 
ffm away from the critical point. As a first step in this direction one can take all loss terms 
due to lightning in the full model into account. In order to have a stochastic process these 
loss terms must be balanced with suitable gain terms. A convenient choice is to attribute 
all gain terms to the completely empty system. Thus, in this generalization of (5.2) we 
introduce transitions from the state | A'^) to the state \L) with weight /(L — N). Then 
the picture (5.3) turns into: 



1 N L-1 L 




(5.8) 



fL 

It is straightforward to obtain the ground state of the Hamiltonian corresponding to (5.8) 
either numerically up to L 10^, or from a continuum approximation. Then one finds for 
the cluster-size distribution 

n{s) ~ (5 9) 

for s ^ 1. It is tempting to interpret this result as a continuously varying critical exponent 
but then one has to introduce a rescaled rate / = fL in order to keep the exponent finite 
in the thermodynamic limit. In any case, the loss processes by lightning strokes in (5.8) 
do not account for the finite correlation length This is surprising because consideration 
of the same processes (with some additional approximations) lead to the correct critical 
behaviour (2.8) of the correlation length [13]- A solution to this puzzle might be the 
following: Once one assumes that there is an upper cutoff, the loss processes by lightning 
strokes fix its critical exponent although these processes alone do not explain the presence 
of this cutoff. 

^) The logarithmic dependence of the relaxation time on the system size is reminiscent 
of results obtained in [25] for certain two-dimensional kinetic spin models with cluster 
dynamics. 
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Other possible modifications of the simphfied Hamiltonian Hg are: 

1) One can introduce death of trees irrespective of their neighbourhood with a rate q, 
thereby modeUng the ignition processes which dominate in the full model if the number 
of trees is small. For / = one then has a simple kinetic model where the spins flip 
independently and n(s) varies exponentially with s. However, this does not give an 
upper cutoff in the general case since then n{s) is well approximated by a sum of the 
results for / = and q = 0, respectively. 

2) One can distribute the loss terms in (5.8) equally over the configurations with fewer 
trees. This amounts to introducing a transition from a state | N^) to a state | Nf) for 
Nf > Ni with rate /. Computing the ground state numerically as before, one finds 
that this model behaves very similarly to that in (5.8). In particular, its cluster-size 
distribution also has the behaviour (5.9). 

Since none of the aforementioned models is able to describe the off-critical behaviour 
correctly, let us finally discuss the accurate projection of the full ffm onto the symmetric 
subspace. This projection can be performed as follows: Consider the matrix element of 
V that leads from A''^ empty places to Nf empty places, i.e. a cluster of Nf — Ni trees is 
struck by lightning. Thus, there must have been L — Nf trees outside the cluster in the 
initial state. This cluster occupies 2 + Nf — N^ places including the two empty places at 
its boundary. This leads to a combinatorial factor {^"^l^j^^^^) in the transition matrix 
element. Translational invar iance yields another factor L and the weight for the process is 
Nf — Ni. Writing this in the normalization (5.1) leads to the matrix element 

{Nf I V I AT.) = ^ ,\ "-"^ ^ (5.10) 

for Nf > Ni. The loss term due to lightning equals {N\V\N) ^ -{L - N) and 
{Nf\V\Ni) = for Nf < Ni. The tree growth has already been correctly accounted 
for in (5.2). Putting this together, one obtains the correct projection of the full Hamilto- 
nian (3.6) onto the symmetric subspace. The corresponding matrix is almost triangular. 

As for the other simplified models, the ground state and thus the cluster-size distribution 
of this Hamiltonian can easily be obtained from a numerical diagonalization for rings with 
up to a few thousands sites. For small / one finds again n{s) ~ for s 1 and 

no crossover to an exponential decay. For / large, n{s) decays rapidly already for small 
s. In this case, the L-dependence of n{s) computed in the symmetric subspace for small 
s becomes small and one can even observe qualitative agreement with simulations of the 
full model. However, even in the case (5.10) the critical behaviour (2.8) for / ^ is 
not reproduced. In summary, the symmetric space seems unable to describe the critical 
behaviour as a function of /, while it does describe the critical correlation functions (3.13) 
very nicely. 
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6. Discussion 



In this paper we have studied the one-dimensional forest-fire model and exhibited several 
new non- integral critical exponents. It remains to be seen to what extent our results are 
relevant for higher dimensions as well. Firstly, it should be checked to a better accuracy 
whether v = vt [14] is indeed true in two dimensions. Another question to be addressed is 
if in higher dimensions the spatial power laws associated with clusters can also be obtained 
by looking at global quantities only. In particular, it remains an open problem if models 
in a symmetric subspace can also describe the critical cluster-size distribution in higher 
dimensions. The two-dimensional version of the ffm introduced in [26] seems to point in 
this direction. In this case, the dynamics is also controlled by a global quantity, namely 
the mean density of trees. 

The fact that the critical density in one dimension is p = 1 leads to the unusual situation 
that all correlation functions at the critical point become independent of the spatial vari- 
ables in the thermodynamic limit. At / = 0, the density in the stationary state equals one 
for all dimensions. Now, a necessary requirement for perturbation theory around / = 
to be of any use is that the density p must be continuous for / — > 0. In one dimension, 
this condition is satisfied, which permitted us to employ perturbation theory and thus led 
to the simplified model. The situation is different in higher dimensions where the critical 
density is smaller than one. Then correlation functions with non-trivial spatial behaviour 
at the critical point are possible, but perturbation theory around / = cannot be used. 

Even if at present our results seem to be specific for one dimension, our simple model in 
the symmetric space shows that spatial power laws can arise from purely global effects. 
A similar observation for certain models of evolution was made recently in [17] where the 
impact of replacing local neighbourhoods by random neighbourhoods was examined. In 
our case, this global dynamics arises naturally at the critical point and is not put in by 
hand. We have also shown that this simplified model has a vanishing gap, i.e. physical 
quantities will relax with a temporal power law to the stationary state. This demonstrates 
that neither power laws in static quantities related to avalanche type processes nor slow 
relaxational modes are a sufficient criterion for criticality in a strict sense, so that the usual 
multi-point correlations should always be studied in addition. 
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Appendix A. A combinatorial identity 

We want to show that 

L-s fL — s — l\ J- 

5(L,.):=5:iS/ = ^ (A.1) 

JV=1 \N-l) 

for < s < L — 1. The proof uses induction on L. First one observes that 

S{L,L- 1)^1 (A.2) 

which gives the start of the induction. Now one can proceed with the following induction 
step 

+ 1, .) = 'f: ' n "^^1:!""" = i + .) . (a.3) 

N=l u=0 

Observing that the r.h.s. of (A.l) satisfies the recurrence relation (A.3) completes the proof 
of this combinatorial identity. 

Appendix B. Perturbation expansion for excitations 

In this appendix we indicate how to compute the excitation spectrum of the ffm using 
quantum-mechanical perturbation theory. First we note that Hq in (3.6) can easily be 
diagonalized. In matrix form the part at one site is given by 

{Ho).= (l . (B.l) 



The right eigenvectors of this matrix are 



(B.2) 



with eigenvalues and 1, respectively. The corresponding left eigenvectors are 

(11), (01). (B.3) 

Now one can treat excited states with standard quantum-mechanical perturbation theory 
where one just has to keep in mind that the left eigenvectors (B.3) are not just the trans- 
posed ones of the right eigenvectors (B.2). From (B.2) and (B.3) it follows immediately 
that the right vectors 

|a;) = -|l...l) + I1...101...1) (B.4) 

X 

and the left vectors 

1 

{y\= Yl (M2...O...I (B.5) 

ii,i2...=0 
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are eigenvectors of Hq with eigenvalue 1. The normahzation in (B.4) and (B.5) is aheady 
{y\x) = 5y^x- Now one easily checks that 



{y\V\x) = 



L-1 



L-1 

















1 











(B.6) 



yx 



The eigenvalues of the matrix (B.6) are straightforward to compute. The result is that H 
has one eigenvalue 1 + f{2L — 1) + 0{f^) and (L — 1) degenerate eigenvalues 1 + f{L — 1) + 
0{p). The latter correspond to the result (3.15) for P ^ 0, and in this case the first-order 
expansion happens to be exact for all / > 0. The new result of the perturbation expansion 
is the eigenvalue 1 + f{2L — 1) + C(/^) which is located in the sector P = 0. 

Similarly, one finds that the L{L — l)/2 vectors 

\xi,X2) = 1 1 ... 1) - I . . . 1 1 . . .) - I . . . 1 1 . . .) + I ... 1 1 ... 1 1 .. .) (B.7) 

Xl X2 Xl X-2 



with Xl < X2 are right eigenvectors of i^o with eigenvalue 2. The corresponding left 
eigenvectors are given by 



{yiM = V (ziZ2...0...0...| (B.8) 

^ ^ ?/l I/O 



. . „ 2/12/2 

ll,l2...=0 



with j/i < 7/2. These vectors are properly normalized, i.e. {yuy^ \ x\-,X'2) = Sy^^xi^y2,x2- ^ 
straightforward computation yields 



{yi,y2\V\xi,X2) = {L-2){l+6y^,a:iSy2,X2)- < 



' {x2 - Xl - 1) if xi<yi,y2< X2 

{yi<xi,x2< y2 
or X2 < yi 
or y2 < Xl 

^ otherwise. 

(B.9) 

It is easy to check that the L vectors \x,x + 1) are eigenvectors of (B.9) with eigenvalue 
(L — 2) . In particular, there is an eigenvalue 2 + f{L — 2) + 0{ f'^) of H for each momentum 
(including the translationally invariant sector) which is compatible with the result (3.15) 
for P^O. 

All first-order corrections discussed so far are positive, but the matrix (B.9) also has 
negative eigenvalues. Since the eigenvalues of a stochastic matrix always have positive 
real parts, this automatically implies that such a negative correction can only be relevant 
for small /. More precisely, one finds eigenvalues 2 — fh{L, P) + 0(/^) where h{L, P) is 
positive around P = tt, symmetric in P — tt and increases for small \P — 7r|. This does 
however not affect the translationally invariant sector, nor is it relevant for the large L 
limit. Nevertheless, the complicated P-dependence of these modes demonstrates that the 
structure of the full relaxation spectrum is highly non-trivial. 
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Appendix C. Continuum approximation 

The model described by Hg in section 5 can be treated in the foUowing continuum approx- 
imation which gives the imaginary part of the spectrum and the eigenfunctions in a simple 
way. 

Denote the eigenvalues by A and the components of the corresponding eigenvector by go, 
gL- Then the eigenvalue problem for (5.2) is given by the following equations: 

fLgo-gi = Xgo: (C.la) 
NgN-{N + l)gN+i = XgN, 1<N<L, (C.lb) 
-fLgo + LgL ^ XgL- (C.lc) 

The solution of the inner equations —{N + l)(^jv+i — gn) — gn = A^jv varies only slowly 
for large N so that we can approximate them by 

-x^g{x) = {X+l)g{x) (C.2) 

for 1 < X < L. Keeping only the leading term in L of the boundary conditions (C.la), 
(C.lc) and eliminating go one finds 

9{L) = ^ . (C.3) 

For the stationary state these equations lead to g{x) = 1/x which corresponds exactly to 
the discrete result gN ^ 1/N. 

The general solution of the difi^erential equation (C.2) is given by 

gix) = a;-(^+i) (C.4) 
and from the boundary condition (C.3) one finds the purely imaginary spectrum 

Afe = ^, keZ. (C.5) 

ln(L) 

The eigenfunctions are plane waves in the variable ln(a;)/ln(L) multiplied by 1/x and 
compare well with those of the discrete problem (C.l), except for small x. 
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